Rare species patterns and results excluding rare species for no-time models
Classifying rare species
Species classification was done for each forest plot and time interval separately. We used the living and dead individuals of the final census from the interval. We used the FuzzyQ clustering algorithm in the ‘FuzzyQ’ R package (Balbuena et al., 2021) to estimate the probability of each species to be common or rare based on species abundance and occupancy in 50x50 m quadrats.
To describe and compare rarity patterns in all forest, we estimated the number of species, number of individuals, and density of rare and common species. For the forests plots with more than 1 census interval, we averaged the estimates across census intervals.
Here, we provide the list of species code and the classification in a Rdata file.
load(here("data/rare_species_classification.Rdata"))
# plot variables including
plots_info <- read.table(here("data/plot_sites_information.csv"), sep=",", header=T)
classi <-classi %>% left_join(plots_info[,c(1,6,9,12,17:23,28)], by=c("fplot" = "ID"))
classi$fplot.time = paste0(classi$fplot,"_", classi$time)rich = classi %>% count(fplot.time,rarity) %>% rename(rich=n)
ntree = classi %>% group_by(fplot.time,rarity) %>%
summarise(ntrees = sum(n))
dat <- rich %>% left_join(ntree, by=c("fplot.time", "rarity")) %>%
pivot_wider(names_from = rarity, values_from =c(rich, ntrees)) %>%
mutate(p.rich_common = rich_common*100/(rich_common+rich_rare),
p.rich_rare = rich_rare*100/(rich_common+rich_rare),
p.ntrees_common = ntrees_common*100/(ntrees_common+ntrees_rare),
p.ntrees_rare = ntrees_rare*100/(ntrees_common+ntrees_rare))
#Mean per forest plots with more than 1 census interval
dat.m <- dat %>% separate(fplot.time, c("fplot", "time")) %>%
group_by(fplot) %>%
summarise_at(vars(rich_common:p.ntrees_rare),mean) %>%
group_by(fplot) %>%
summarise_at(vars(rich_common:p.ntrees_rare),round,digits=0) %>%
select(fplot,rich_common, p.rich_common, rich_rare, p.rich_rare,
ntrees_common, p.ntrees_common, ntrees_rare, p.ntrees_rare)dat.m %>%
htmlTable(header=c("Forest" ,rep(c("N", "%"), 4)),
cgroup = list(c("", "Species richness", "Number of trees"),
c("", "Common", "Rare", "Common", "rare")),
n.cgroup = list(c(1,4,4),
c(1,2,2,2,2)),caption="Table S4.1. Number and percentages of species and trees classified as common or rare per forest plot. See Table S1.1 and S1.2 for forest plot names and information. For forest plots with more than 1 census interval, we average the values across intervals."
)| Table S4.1. Number and percentages of species and trees classified as common or rare per forest plot. See Table S1.1 and S1.2 for forest plot names and information. For forest plots with more than 1 census interval, we average the values across intervals. | |||||||||||||
| Species richness | Number of trees | ||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Common | Rare | Common | rare | ||||||||||
| Forest | N | % | N | % | N | % | N | % | |||||
| 1 | ama | 456 | 35 | 839 | 65 | 98878 | 90 | 11578 | 10 | ||||
| 2 | bci | 135 | 42 | 185 | 58 | 334275 | 96 | 12477 | 4 | ||||
| 3 | edo | 122 | 29 | 298 | 71 | 168253 | 96 | 7004 | 4 | ||||
| 4 | fus | 55 | 50 | 55 | 50 | 143813 | 98 | 2376 | 2 | ||||
| 5 | idc | 74 | 53 | 65 | 47 | 34651 | 95 | 1675 | 5 | ||||
| 6 | kor | 186 | 40 | 282 | 60 | 335392 | 94 | 22608 | 6 | ||||
| 7 | lam | 539 | 38 | 869 | 62 | 376389 | 86 | 61135 | 14 | ||||
| 8 | ldw | 14 | 38 | 23 | 62 | 28004 | 96 | 1257 | 4 | ||||
| 9 | len | 108 | 29 | 262 | 71 | 145467 | 97 | 4296 | 3 | ||||
| 10 | lpl | 108 | 45 | 132 | 55 | 120516 | 96 | 4735 | 4 | ||||
| 11 | luq | 55 | 37 | 95 | 63 | 69233 | 96 | 2769 | 4 | ||||
| 12 | mos | 96 | 34 | 182 | 66 | 157032 | 95 | 7652 | 5 | ||||
| 13 | pas | 375 | 44 | 478 | 56 | 356955 | 92 | 31286 | 8 | ||||
| 14 | scbi | 25 | 35 | 47 | 65 | 38366 | 97 | 1248 | 3 | ||||
| 15 | serc | 20 | 27 | 54 | 73 | 27031 | 96 | 1204 | 4 | ||||
| 16 | sin | 118 | 50 | 118 | 50 | 207356 | 94 | 12196 | 6 | ||||
| 17 | ucsc | 10 | 32 | 21 | 68 | 8738 | 96 | 394 | 4 | ||||
| 18 | wab | 11 | 28 | 28 | 72 | 49562 | 93 | 3558 | 7 | ||||
| 19 | wfdp | 8 | 31 | 18 | 69 | 28948 | 97 | 831 | 3 | ||||
| 20 | wyw | 5 | 20 | 20 | 80 | 19300 | 96 | 732 | 4 | ||||
| 21 | zof | 3 | 23 | 10 | 77 | 75368 | 99 | 381 | 1 | ||||
classi %>% group_by(rarity) %>%
mutate(dens.ha = n/size_data) %>%
summarise(minD = min(dens.ha),
maxD = max(dens.ha),
meanD = mean(dens.ha),
sdD = sd(dens.ha),
medianD = median(dens.ha),
q90D = quantile(dens.ha, 0.9),
q95D = quantile(dens.ha, 0.95)) %>%
kable(digits=2, caption = "Table S4.2. Summary values for the density of trees (N/ha) classified as common or rare for the 21 forest plots. SD - standard deviation, Quant - quantiles.")| rarity | minD | maxD | meanD | sdD | medianD | q90D | q95D |
|---|---|---|---|---|---|---|---|
| common | 1.32 | 3840.80 | 35.24 | 128.01 | 11.26 | 63.74 | 109.80 |
| rare | 0.02 | 297.32 | 1.39 | 5.14 | 0.55 | 2.84 | 4.39 |
dat.m2 <- dat.m %>%
left_join(plots_info[,c(1,6,9,12,17:23,28)], by=c("fplot" = "ID")) %>%
mutate(tot.rich = rich_common + rich_rare,
tot.n = ntrees_common+ntrees_rare)prich <- dat.m2 %>%
pivot_longer(c(p.rich_common, p.rich_rare),names_to = "rarity",
values_to = "richness") %>%
ggplot(aes(x=fct_reorder(fplot, tot.rich, max),
y=richness, fill=rarity)) +
geom_col() +
xlab("") + ylab("Species richness (%)") +
labs(tag="A")+
scale_fill_discrete("Species",
labels=c("common", "rare"))+
theme(axis.text.x = element_text(angle = 45, hjust=1))
p.ntree <- dat.m2 %>%
pivot_longer(c(p.ntrees_common, p.ntrees_rare),names_to = "rarity",
values_to = "ntrees") %>%
ggplot(aes(x=fct_reorder(fplot, Nsp_literature, max),
y=ntrees, fill=rarity)) +
geom_col()+
labs(tag="B")+
xlab("") + ylab("Number of trees (%)") +
scale_fill_discrete("Species",
labels=c("common", "rare"))+
theme(axis.text.x = element_text(angle = 45, hjust=1),
legend.position = "none")
prich/p.ntree + plot_layout(guides = "collect")Figure S4.1. Percentages of rare and common species and number of trees per forest plot. Forest plots are arranged by absolute latitude (from the largest to the smallest). See Table S1.2 for plots names.
library(wesanderson) #better colour palette
# Gradient color for latitud
pal <- wes_palette("Zissou1",21, type = "continuous")[21:1] # azul frio-temperadodat.m2 %>%
ggplot(aes(x=tot.rich, y=p.rich_rare))+
scale_x_log10()+
geom_point(aes( col=abs(lat)), size=3) +
scale_color_gradientn(name="Latitude \n (abs)",colors=pal) +
geom_smooth(method=mgcv::gam,se=T, formula = y ~ s(x, bs = "cs")) +
xlab("Number of species") +
ylab("Percentage of rare species") Figure S4.2. Percentage of rare species against the number of species (log10 scale) for 21 worldwide distributed forest plots. Blue line and grey area are the fitted results and confidence intervals for a generalised additive model showing a decrease in the percentage of rare species with the number of species but only for forests with less than 100 species. Model’s adjusted R2 = 0.30. Each forest plot is coloured by the latitude in absolute values.
mod <- mgcv::gam(p.rich_rare ~ s(log(tot.rich), bs="cs"), data=dat.m2, method = "REML")
summary(mod)##
## Family: gaussian
## Link function: identity
##
## Formula:
## p.rich_rare ~ s(log(tot.rich), bs = "cs")
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 63.810 1.632 39.09 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(log(tot.rich)) 1.803 9 0.958 0.0168 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.301 Deviance explained = 36.4%
## -REML = 72.24 Scale est. = 55.953 n = 21
No time models
Running models with different datasets:
all: original data, all species all individuals
exclude: exclude rare species
regroup: regrouping observations of rare species in ONE species group
Using only the 5x5 m quadrat scale
#Figures colors/labels
cores <- c("#c00000", "#b666d2", "#70ad47", "#ffc000")
labvit <- as_labeller(c(grow = "Growth", mort="Mortality",
rec="Recruitment"))
grplab = as_labeller(c(quadrat = "Space", `quadrat:sp`= "Space X Species",
sp="Species", Residual = "residual"))
grpvit <- as_labeller(c(grow = "Growth", mort="Mortality",
rec="Recruitment",
quadrat = "Space", `quadrat:sp`= "Species X Space",
sp="Species", Residual = "Residual"))Data
Growth
grp <- c("all","exclude", "regroup")
tudo <- list()
path1 = here("models_outputs/models_notime/")
for (i in 1:length(grp)){
path <- paste0(path1, grp[i], "/grow/table")
files = list.files(path)
all <- lapply(files, function(x) get(assign(x, load(paste0(path, "/", x)))))
names(all) <- substr(files,1,nchar(files)-6)
tudo[[i]] <- all
}
names(tudo) <- grp
tes <- map(tudo,bind_rows, .id="id") %>% bind_rows() %>% filter(q_size=="quad_5")
#tes %>% count(data,fplot) %>% mutate(n=n/8)
table(tes$fplot,tes$data)/5##
## all exclude regroup
## ama 1 1 1
## bci 6 6 6
## edo 2 2 2
## fus 3 3 3
## idc 1 1 1
## kor 1 1 1
## lam 3 3 3
## ldw 1 1 1
## len 2 2 2
## lpl 1 1 1
## luq 3 3 3
## mos 2 2 2
## pas 4 4 4
## scbi 1 1 1
## serc 1 1 1
## sin 2 2 2
## ucsc 2 2 2
## wab 2 2 2
## wfdp 1 1 1
## wyw 3 3 3
## zof 1 1 1
#remove intercept to include again as a separate column
intercepts <- tes %>% filter(term == "intercept")
names(intercepts)[9] = "intercept"
grow <- tes %>% filter(term != "intercept") %>%
left_join(intercepts[,c(1:4,9)], by= c("id", "data", "fplot", "time")) %>%
group_by(id, data, fplot, time) %>%
mutate(VPC = variance/sum(variance)) %>% ungroup() %>%
mutate(term = fct_relevel(term, "quadrat", "quadrat:sp", "sp", "Residual")) %>%
as.data.frame()
# divide the sd of the terms by the mean growth (intercept) for each dataset.
grow$stdev = grow$Estimate/grow$interceptMort
grp <- c("all","exclude", "regroup")
tudo <- list()
path1 = here("models_outputs/models_notime/")
for (i in 1:length(grp)){
path <- paste0(path1,grp[i], "/mort/table")
files = list.files(path)
all <- lapply(files, function(x) get(assign(x, load(paste0(path, "/", x)))))
names(all) <- substr(files,1,nchar(files)-6)
tudo[[i]] <- all
}
names(tudo) <- grp
tes <- map(tudo,bind_rows, .id="id") %>% bind_rows() %>% filter(q_size=="quad_5")
#tes %>% count(data,fplot) %>% mutate(n=n/8)
table(tes$fplot,tes$data)/5##
## all exclude regroup
## ama 1 1 1
## bci 7 7 7
## edo 2 2 2
## fus 3 3 3
## idc 1 1 1
## kor 1 1 1
## lam 3 3 3
## ldw 1 1 1
## len 2 2 2
## lpl 1 1 1
## luq 3 3 3
## mos 2 2 2
## pas 5 5 5
## scbi 1 1 1
## serc 1 1 1
## sin 2 2 2
## ucsc 2 2 2
## wab 2 2 2
## wfdp 1 1 1
## wyw 3 3 3
## zof 1 1 1
#remove intercept to include again as a separate column
intercepts <- tes %>% filter(term == "intercept")
names(intercepts)[9] = "intercept"
mort <- tes %>% filter(term != "intercept") %>%
left_join(intercepts[,c(1:4,9)], by= c("id", "data", "fplot", "time")) %>%
group_by(id, data, fplot, time) %>%
mutate(VPC = variance/sum(variance)) %>% ungroup() %>%
mutate(term = fct_relevel(term, "quadrat", "quadrat:sp", "sp", "Residual")) %>%
as.data.frame()
mort$stdev = mort$EstimateRec
grp <- c("all","exclude", "regroup")
tudo <- list()
path1 = here("models_outputs/models_notime/")
for (i in 1:length(grp)){
path <- paste0(path1,grp[i], "/rec/table")
files = list.files(path)
all <- lapply(files, function(x) get(assign(x, load(paste0(path, "/", x)))))
names(all) <- substr(files,1,nchar(files)-6)
tudo[[i]] <- all
}
names(tudo) <- grp
tes <- map(tudo,bind_rows, .id="id") %>% bind_rows() %>% filter(q_size=="quad_5")
#tes %>% count(data,fplot) %>% mutate(n=n/8)
table(tes$fplot,tes$data)/5##
## all exclude regroup
## bci 7 7 7
## edo 2 2 2
## fus 3 3 3
## idc 1 1 1
## kor 1 1 1
## lam 3 3 3
## ldw 1 1 1
## len 2 2 2
## lpl 1 1 1
## luq 3 3 3
## mos 2 2 2
## pas 5 5 5
## scbi 1 1 1
## serc 1 1 1
## sin 2 2 2
## ucsc 2 2 2
## wab 2 2 2
## wfdp 1 1 1
## wyw 3 3 3
## zof 1 1 1
#remove intercept to include again as a separate column
intercepts <- tes %>% filter(term == "intercept")
names(intercepts)[9] = "intercept"
rec <- tes %>% filter(term != "intercept") %>%
left_join(intercepts[,c(1:4,9)], by= c("id", "data", "fplot", "time")) %>%
group_by(id, data, fplot, time) %>%
mutate(VPC = variance/sum(variance)) %>% ungroup() %>%
mutate(term = fct_relevel(term, "quadrat", "quadrat:sp", "sp", "Residual")) %>%
as.data.frame()
rec$stdev = rec$EstimateExcluding recruitment for Wytham Woods, as in “all data”.
Combining results.
Comparing groups
comb %>%
ggplot(aes(x=data, y=VPC)) +
geom_col(aes(fill=term, color=term), alpha=0.7) +
scale_fill_manual(values=cores) +
scale_color_manual(values=cores) +
facet_grid(vital~fplot+time) +
theme(axis.text.x = element_text(angle=45, hjust=1))mean over fplots
comb %>% group_by(vital,data,term) %>% summarise(VPC=mean(VPC)) %>%
ggplot(aes(x=data, y=VPC)) +geom_col(aes(fill=term, color=term), alpha=0.7) +
scale_fill_manual(values=cores) +
scale_color_manual(values=cores) +
facet_grid(~vital)+
theme(axis.text.x = element_text(angle=45, hjust=1))+
labs(tag="a)")comb %>%
ggplot(aes(x=data,y=VPC, col=data)) + geom_boxplot()+
facet_grid(vital~term) +
theme(axis.text.x = element_text(angle=45,hjust=1),
legend.position = "none")comb %>% mutate(xis = as.numeric(as.factor(data))) %>%
ggplot(aes(x=xis,y=stdev, col=paste(fplot,time))) + geom_line() +
facet_grid(vital~term) +
scale_x_continuous(breaks = 1:3, labels=levels(comb$data),
name="") +
theme(axis.text.x = element_text(angle=45, hjust=1),
legend.position = "none")comb %>% mutate(xis = as.numeric(as.factor(data))) %>%
ggplot(aes(x=xis,y=VPC, col=paste(fplot,time))) + geom_line() +
facet_grid(vital~term) +
scale_x_continuous(breaks = 1:3, labels=levels(comb$data),
name="") +
theme(axis.text.x = element_text(angle=45, hjust=1),
legend.position = "none")Difference to all data
Absolute difference
dif <- comb %>%
pivot_wider(id_cols = c(vital,fplot,time,term),
names_from = data,
values_from = c(stdev,VPC, richness)) %>%
mutate(stdev.dif_exclude = stdev_exclude - stdev_all,
VPC.dif_exclude = VPC_exclude - VPC_all,
stdev.dif_regroup = stdev_regroup - stdev_all,
VPC.dif_regroup = VPC_regroup - VPC_all) %>%
select(vital,fplot,time,term,richness_exclude,
stdev.dif_exclude,
stdev.dif_regroup, VPC.dif_exclude, VPC.dif_regroup) %>%
pivot_longer(6:9, names_to = c("index", "data"), names_sep = "_",
#names_pattern = "(.*_).(.*)",
values_to = "value") %>%
pivot_wider(names_from = index, values_from = value)dif %>%
ggplot(aes(x=data, y=stdev.dif)) +
facet_grid(vital~term, scales= "free_x") +
geom_boxplot() +
geom_quasirandom(aes(col=fplot)) +
ylab("SD no.rare - SD all") + xlab("") +
geom_hline(yintercept = 0, linetype= "dotted") +
ggtitle("Difference in SD") +
theme(legend.position = "right",
panel.background = element_rect(colour="lightgray"),
axis.text.x = element_text(angle=45, hjust=1))dif %>%
ggplot(aes(x=data, y=VPC.dif)) +
facet_grid(vital~term, scales= "free_x") +
geom_boxplot() +
geom_quasirandom(aes(col=fplot)) +
ylab("VPC no.rare - VPC all") + xlab("") +
geom_hline(yintercept = 0, linetype= "dotted") +
ggtitle("Difference in VPC") +
theme(legend.position = "right",
panel.background = element_rect(colour="lightgray"),
axis.text.x = element_text(angle=45, hjust=1))Relative difference SD
difr <- comb %>%
pivot_wider(id_cols = c(vital,fplot,time,term),
names_from = data,
values_from = c(stdev,VPC, richness)) %>%
mutate(stdev.dif_exclude = (stdev_exclude - stdev_all)*100/stdev_all,
stdev.dif_regroup = (stdev_regroup - stdev_all)*100/stdev_all) %>%
select(vital,fplot,time,term,richness_exclude,
stdev.dif_exclude,stdev.dif_regroup) %>%
pivot_longer(6:7, names_to = "data",
values_to = "stdev.dif") %>%
mutate(data = substr(data,11, nchar(data)))difr %>%
ggplot(aes(x=data, y=stdev.dif)) +
facet_grid(vital~term, scales= "free_x") +
geom_boxplot() +
geom_quasirandom(aes(col=fplot)) +
ylab("Proportional difference (%)") + xlab("") +
geom_hline(yintercept = 0, linetype= "dotted") +
ggtitle("Relative difference in SD") +
theme(legend.position = "right",
panel.background = element_rect(colour="lightgray"),
axis.text.x = element_text(angle=45, hjust=1))difr %>% group_by(vital, term, data) %>% summarise(mean=mean(stdev.dif, na.rm=T)) %>%
pivot_wider(names_from = term, values_from = mean)## # A tibble: 6 × 6
## # Groups: vital [3]
## vital data quadrat `quadrat:sp` sp Residual
## <chr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 grow exclude 11.0 8.47 -11.8 13.2
## 2 grow regroup 8.71 3.97 -12.7 13.9
## 3 mort exclude 0.0278 -1.03 -14.9 0
## 4 mort regroup -2.62 -0.591 -16.4 0
## 5 rec exclude -0.565 -0.299 -9.26 0
## 6 rec regroup 0.512 -0.768 -10.4 0
table
Relative differences in SD
difr %>% group_by(vital, term, data) %>%
summarise(mean=mean(stdev.dif, na.rm=T)) %>%
pivot_wider(names_from = c(vital,data), values_from = mean) %>%
kable(digits=1)| term | grow_exclude | grow_regroup | mort_exclude | mort_regroup | rec_exclude | rec_regroup |
|---|---|---|---|---|---|---|
| quadrat | 11.0 | 8.7 | 0.0 | -2.6 | -0.6 | 0.5 |
| quadrat:sp | 8.5 | 4.0 | -1.0 | -0.6 | -0.3 | -0.8 |
| sp | -11.8 | -12.7 | -14.9 | -16.4 | -9.3 | -10.4 |
| Residual | 13.2 | 13.9 | 0.0 | 0.0 | 0.0 | 0.0 |
Mean absolute differences in VPC
dif %>% group_by(vital, term, data) %>%
summarise(mean=mean(VPC.dif, na.rm=T)) %>%
pivot_wider(names_from = c(vital,data), values_from = mean) -> write
kable(write)| term | grow_exclude | grow_regroup | mort_exclude | mort_regroup | rec_exclude | rec_regroup |
|---|---|---|---|---|---|---|
| quadrat | 0.0048918 | 0.0039062 | 0.0081885 | 0.0061325 | 0.0069821 | 0.0130636 |
| quadrat:sp | 0.0079957 | 0.0028262 | 0.0139204 | 0.0162427 | 0.0076814 | 0.0044177 |
| sp | -0.0831334 | -0.0908756 | -0.0704275 | -0.0754252 | -0.0321646 | -0.0368164 |
| Residual | 0.0702459 | 0.0841432 | 0.0483186 | 0.0530499 | 0.0175012 | 0.0193351 |
Comparing with NUMBER of rare species
dif %>%
ggplot(aes(x=richness_exclude, y=stdev.dif, shape=data)) +
geom_point(aes(col=fplot)) +
geom_smooth(method="lm", se=F, aes(linetype=data)) +
facet_grid(vital~term) +
ylab("Difference in SD to original data") +
xlab("log N of rare species") +
geom_hline(yintercept = 0, linetype= "dotted") +
ggtitle("Difference in SD") +
scale_x_log10()+
theme(axis.text.x = element_text(angle=45, hjust=1),
)difr %>%
ggplot(aes(x=richness_exclude, y=stdev.dif, shape=data)) +
geom_point(aes(col=fplot)) +
geom_smooth(method="lm", se=F, aes(linetype=data)) +
facet_grid(vital~term) +
ylab("Relative difference in SD to original data") +
xlab("log N of rare species") +
geom_hline(yintercept = 0, linetype= "dotted") +
ggtitle("Proportional difference in SD to original data") +
theme(axis.text.x = element_text(angle=45, hjust=1),
panel.background = element_rect(fill="white", color="black"))pm <- dif %>%
ggplot(aes(x=richness_exclude, y=VPC.dif, shape=data)) +
geom_point(aes(col=fplot)) +
facet_grid(vital~term) +
geom_smooth(method="lm", se=F, aes(linetype=data)) +
ylab(" Absolute difference in VPC to original data") +
xlab("log N of rare species") +
geom_hline(yintercept = 0, linetype= "dotted") +
theme(axis.text.x = element_text(angle=45, hjust=1),
panel.background = element_rect(fill="white", color="black"))
pmDirichlet regression for exclude results with species richness
Richness by rarefaction at the 6ha plot size (sampling unit 20x20m)
load(here("data", "rarefaction.curves.Rdata"))
rich.rare <- rare[rare$sites ==150,2:3]
colnames(rich.rare)[1] <- "richness.rarefaction"
comb <- comb %>% left_join(rich.rare, by="fplot")Latitude
load(here("data", "plots_structure.Rdata"))
comb <- comb %>% left_join(plots.structure[,1:2], "fplot")#mean for each forest
mcomb <- comb %>% group_by(vital,data, fplot,term, lat) %>%
summarise(VPC = mean(VPC),
stdev = mean(stdev),
richness.rarefaction = mean(richness.rarefaction))GROWTH
predirig <- mcomb %>%
filter(vital == "grow", data == "exclude") %>%
select(fplot, term, VPC, richness.rarefaction,lat)
diridatag <- predirig %>%
pivot_wider(names_from = term, values_from = VPC) %>% ungroup() %>%
mutate(log.rich = log(richness.rarefaction),
log.rich.o=log.rich) %>%
mutate_at(vars( log.rich), scale)
vpcg <- DR_data(diridatag[,c("quadrat", "sp", "quadrat:sp", "Residual")])
#plot(vpc)## Call:
## DirichReg(formula = vpcg ~ log.rich, data = diridatag, model = "alternative",
## base = 4)
##
## Standardized Residuals:
## Min 1Q Median 3Q Max
## quadrat -0.9058 -0.4992 -0.2701 -0.0166 0.9614
## sp -1.5134 -0.7438 -0.3542 0.2168 3.7107
## quadrat:sp -1.7533 -0.5644 -0.2376 0.3694 2.6473
## Residual -2.0096 -0.8659 0.3287 1.1003 1.8630
##
## MEAN MODELS:
## ------------------------------------------------------------------
## Coefficients for variable no. 1: quadrat
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.6205 0.1744 -15.026 < 2e-16 ***
## log.rich 0.4740 0.1675 2.829 0.00466 **
## ------------------------------------------------------------------
## Coefficients for variable no. 2: sp
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.09746 0.09944 -11.036 <2e-16 ***
## log.rich -0.20184 0.10053 -2.008 0.0447 *
## ------------------------------------------------------------------
## Coefficients for variable no. 3: quadrat:sp
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.44533 0.11231 -12.869 <2e-16 ***
## log.rich -0.08056 0.11309 -0.712 0.476
## ------------------------------------------------------------------
## Coefficients for variable no. 4: Residual
## - variable omitted (reference category) -
## ------------------------------------------------------------------
##
## PRECISION MODEL:
## ------------------------------------------------------------------
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.4105 0.1789 19.06 <2e-16 ***
## ------------------------------------------------------------------
## Significance codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Log-likelihood: 102.9 on 7 df (38 BFGS + 1 NR Iterations)
## AIC: -191.7, BIC: -184.4
## Number of Observations: 21
## Links: Logit (Means) and Log (Precision)
## Parametrization: alternative
newdata <- expand.grid(log.rich = seq(min(diridatag$log.rich),
max(diridatag$log.rich), length.out=10))
pred <- predict(grow2, newdata = newdata, se=T)
confint(grow2)##
## 95% Confidence Intervals (original form)
##
## - Beta-Parameters:
## Variable: quadrat
## 2.5% Est. 97.5%
## (Intercept) -2.962 -2.621 -2.279
## log.rich 0.146 0.474 0.802
##
## Variable: sp
## 2.5% Est. 97.5%
## (Intercept) -1.292 -1.097 -0.903
## log.rich -0.399 -0.202 -0.005
##
## Variable: quadrat:sp
## 2.5% Est. 97.5%
## (Intercept) -1.665 -1.445 -1.225
## log.rich -0.302 -0.081 0.141
##
## Variable: Residual
## variable omitted
##
## - Gamma-Parameters
## 2.5% Est. 97.5%
## (Intercept) 3.06 3.41 3.76
colnames(pred) <- colnames(vpcg)
newdata$log.rich.o <- newdata$log.rich*sd(diridatag$log.rich.o) +
mean(diridatag$log.rich.o)
newdata$rich.o <- exp(newdata$log.rich.o)
newpredg <- cbind(newdata,pred) %>% pivot_longer(`quadrat`:`Residual`, names_to="term", values_to="pred")
newpredg %>% mutate(term = fct_relevel(term, "quadrat", "quadrat:sp", "sp",
"Residual")) %>%
ggplot(aes(x=rich.o, y=pred, col=term)) +
geom_line() + facet_grid(~term) +
geom_point(data=predirig, aes(x=richness.rarefaction, y=VPC, col=term)) +
scale_x_log10() +
ggtitle("Mod log.rich")Residuals
resid <- residuals(grow2, type = "standardized")
resid <- data.frame(resid[,1:4])
resid$log.rich <- diridatag$log.rich
resid <- resid %>% pivot_longer(1:4, names_to="term", values_to="resid")
ggplot(resid, aes(x=log.rich, y=resid)) + geom_point() +
facet_grid(~term)MORTALITY
predirim <- mcomb %>%
filter(vital == "mort", data == "exclude") %>%
select(fplot, term, VPC, richness.rarefaction,lat)
diridatam <- predirim %>%
pivot_wider(names_from = term, values_from = VPC) %>% ungroup() %>%
mutate(log.rich = log(richness.rarefaction),
log.rich.o=log.rich) %>%
mutate_at(vars( log.rich), scale)
vpcm <- DR_data(diridatam[,c("quadrat", "sp", "quadrat:sp", "Residual")])
#plot(vpc)## Call:
## DirichReg(formula = vpcm ~ log.rich, data = diridatam, model = "alternative",
## base = 4)
##
## Standardized Residuals:
## Min 1Q Median 3Q Max
## quadrat -1.2846 -0.6061 -0.0454 0.4078 1.8812
## sp -1.9682 -0.8456 0.1570 0.9680 2.6456
## quadrat:sp -1.8027 -0.6593 -0.0102 0.2449 1.6414
## Residual -2.9957 -0.2679 0.0598 0.6229 1.6857
##
## MEAN MODELS:
## ------------------------------------------------------------------
## Coefficients for variable no. 1: quadrat
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.6901 0.1375 -12.295 < 2e-16 ***
## log.rich 0.4103 0.1473 2.786 0.00534 **
## ------------------------------------------------------------------
## Coefficients for variable no. 2: sp
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.87149 0.10008 -8.708 <2e-16 ***
## log.rich -0.14106 0.09865 -1.430 0.153
## ------------------------------------------------------------------
## Coefficients for variable no. 3: quadrat:sp
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.2178 0.1135 -10.73 <2e-16 ***
## log.rich -0.2471 0.1217 -2.03 0.0424 *
## ------------------------------------------------------------------
## Coefficients for variable no. 4: Residual
## - variable omitted (reference category) -
## ------------------------------------------------------------------
##
## PRECISION MODEL:
## ------------------------------------------------------------------
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.3788 0.1756 19.24 <2e-16 ***
## ------------------------------------------------------------------
## Significance codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Log-likelihood: 88.64 on 7 df (43 BFGS + 1 NR Iterations)
## AIC: -163.3, BIC: -156
## Number of Observations: 21
## Links: Logit (Means) and Log (Precision)
## Parametrization: alternative
newdata <- expand.grid(log.rich = seq(min(diridatam$log.rich),
max(diridatam$log.rich), length.out=10))
pred <- predict(mort2, newdata = newdata, se=T)
confint(mort2)##
## 95% Confidence Intervals (original form)
##
## - Beta-Parameters:
## Variable: quadrat
## 2.5% Est. 97.5%
## (Intercept) -1.960 -1.69 -1.421
## log.rich 0.122 0.41 0.699
##
## Variable: sp
## 2.5% Est. 97.5%
## (Intercept) -1.068 -0.871 -0.675
## log.rich -0.334 -0.141 0.052
##
## Variable: quadrat:sp
## 2.5% Est. 97.5%
## (Intercept) -1.440 -1.218 -0.995
## log.rich -0.486 -0.247 -0.009
##
## Variable: Residual
## variable omitted
##
## - Gamma-Parameters
## 2.5% Est. 97.5%
## (Intercept) 3.04 3.38 3.72
colnames(pred) <- colnames(vpcm)
newdata$log.rich.o <- newdata$log.rich*sd(diridatam$log.rich.o) +
mean(diridatam$log.rich.o)
newdata$rich.o <- exp(newdata$log.rich.o)
newpredm <- cbind(newdata,pred) %>% pivot_longer(`quadrat`:`Residual`, names_to="term", values_to="pred")
newpredm %>% mutate(term = fct_relevel(term, "quadrat", "quadrat:sp", "sp",
"Residual")) %>%
ggplot(aes(x=rich.o, y=pred, col=term)) +
geom_line() + facet_grid(~term) +
geom_point(data=predirim, aes(x=richness.rarefaction, y=VPC, col=term)) +
scale_x_log10() +
ggtitle("Mod log.rich")Residuals
resid <- residuals(mort2, type = "standardized")
resid <- data.frame(resid[,1:4])
resid$log.rich <- diridatam$log.rich
resid <- resid %>% pivot_longer(1:4, names_to="term", values_to="resid")
ggplot(resid, aes(x=log.rich, y=resid)) + geom_point() +
facet_grid(~term)RECRUITMENT
predirir <- mcomb %>%
filter(vital == "rec", data == "exclude") %>%
select(fplot, term, VPC, richness.rarefaction,lat)
diridatar <- predirir %>%
pivot_wider(names_from = term, values_from = VPC) %>% ungroup() %>%
mutate(log.rich = log(richness.rarefaction),
log.rich.o=log.rich) %>%
mutate_at(vars( log.rich), scale)
vpcr <- DR_data(diridatar[,c("quadrat", "sp", "quadrat:sp", "Residual")])
#plot(vpc)## Call:
## DirichReg(formula = vpcr ~ log.rich, data = diridatar, model = "alternative",
## base = 4)
##
## Standardized Residuals:
## Min 1Q Median 3Q Max
## quadrat -1.5124 -0.6320 -0.2358 0.8198 2.1508
## sp -2.8666 -0.2882 0.1438 0.7867 1.6337
## quadrat:sp -1.4254 -0.6122 -0.2100 0.3186 2.8831
## Residual -1.5080 -0.6407 -0.0523 0.5100 3.1244
##
## MEAN MODELS:
## ------------------------------------------------------------------
## Coefficients for variable no. 1: quadrat
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.586194 0.129600 -4.523 6.09e-06 ***
## log.rich 0.005999 0.140757 0.043 0.966
## ------------------------------------------------------------------
## Coefficients for variable no. 2: sp
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.08768 0.11441 -0.766 0.443
## log.rich -0.71600 0.11995 -5.969 2.38e-09 ***
## ------------------------------------------------------------------
## Coefficients for variable no. 3: quadrat:sp
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.70344 0.13331 -5.277 1.32e-07 ***
## log.rich -0.07571 0.14027 -0.540 0.589
## ------------------------------------------------------------------
## Coefficients for variable no. 4: Residual
## - variable omitted (reference category) -
## ------------------------------------------------------------------
##
## PRECISION MODEL:
## ------------------------------------------------------------------
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.2418 0.1825 17.77 <2e-16 ***
## ------------------------------------------------------------------
## Significance codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Log-likelihood: 71.27 on 7 df (43 BFGS + 1 NR Iterations)
## AIC: -128.5, BIC: -121.9
## Number of Observations: 19
## Links: Logit (Means) and Log (Precision)
## Parametrization: alternative
newdata <- expand.grid(log.rich = seq(min(diridatar$log.rich),
max(diridatar$log.rich), length.out=10))
pred <- predict(rec2, newdata = newdata, se=T)
colnames(pred) <- colnames(vpcr)
newdata$log.rich.o <- newdata$log.rich*sd(diridatar$log.rich.o) +
mean(diridatar$log.rich.o)
newdata$rich.o <- exp(newdata$log.rich.o)
newpredr <- cbind(newdata,pred) %>% pivot_longer(`quadrat`:`Residual`, names_to="term", values_to="pred")
newpredr %>% mutate(term = fct_relevel(term, "quadrat", "quadrat:sp", "sp",
"Residual")) %>%
ggplot(aes(x=rich.o, y=pred, col=term)) +
geom_line() + facet_grid(~term) +
geom_point(data=predirir, aes(x=richness.rarefaction, y=VPC, col=term)) +
scale_x_log10() +
ggtitle("Mod log.rich")Residuals
resid <- residuals(rec2, type = "standardized")
resid <- data.frame(resid[,1:4])
resid$log.rich <- diridatar$log.rich
resid <- resid %>% pivot_longer(1:4, names_to="term", values_to="resid")
ggplot(resid, aes(x=log.rich, y=resid)) + geom_point() +
facet_grid(~term)FIGURE
Calculated prediction interval
source(here("scripts/prediction_intervals_dirichlet_exclude.R"), local=T)
#load("results/prediction_intervals_dirichlet_exclude.Rdata") # quantspredis <- bind_rows(list(grow=newpredg, mort=newpredm, rec=newpredr),
.id="vital") %>%
mutate(term = fct_relevel(term, "quadrat", "quadrat:sp", "sp","Residual"))
pontos <- bind_rows(list(grow=predirig, mort= predirim, rec=predirir), .id="vital")%>%
mutate(term = fct_relevel(term, "quadrat", "quadrat:sp", "sp", "Residual"))
quants <- bind_rows(list(grow=quantg,mort=quantm, rec=quantr), .id="vital")%>%
mutate(term = fct_relevel(term, "quadrat", "quadrat:sp", "sp", "Residual"))# pvalues
res <- data.frame(vital = c("grow", "mort", "rec"),
P = NA,
term = "Residual")
test <- summary(grow2)
vals <- as.data.frame(test$coef.mat)[1:6,]
vals$coef <- names(test$coefficients)[1:6]
valsg <- vals[ grep("log.rich", vals$coef),] %>%
mutate(term=c("quadrat", "sp", "quadrat:sp"))
test <- summary(mort2)
vals <- as.data.frame(test$coef.mat)[1:6,]
vals$coef <- names(test$coefficients)[1:6]
valsm <- vals[ grep("log.rich", vals$coef),] %>%
mutate(term=c("quadrat", "sp", "quadrat:sp"))
test <- summary(rec2)
vals <- as.data.frame(test$coef.mat)[1:6,]
vals$coef <- names(test$coefficients)[1:6]
valsr <- vals[ grep("log.rich", vals$coef),] %>%
mutate(term=c("quadrat", "sp", "quadrat:sp"))
pvals <- bind_rows(list(grow=valsg,mort=valsm,rec=valsr), .id="vital") %>%
dplyr::select(vital,`Pr(>|z|)`, term) %>% rename(P = `Pr(>|z|)`)
pvals <- rbind(pvals,res) %>%
mutate(term = fct_relevel(term, "quadrat", "quadrat:sp", "sp", "Residual"))
pvals$x <- 300
pvals$y <- 0.70
pvals$sig <- paste0("p = ", round(pvals$P,3))
#pvals$sig[pvals$P >0.0166] <- "ns"
pvals$sig[pvals$term=="Residual"] <- ""
pvals$sig[pvals$sig=="p = 0"] <- "p < 0.0001"library(wesanderson) #better colour palette
# Gradient color para latitutde
pal <- wes_palette("Zissou1",20, type = "continuous")[20:1] # azul frio-temperadofdiri_lat <-ggplot(predis, aes(x=rich.o, y=pred)) +
geom_line(size=1) +
facet_grid(vital~term, labeller=grpvit) +
geom_smooth(data=quants,aes(x=rich.o, y=lower), se=F, size=0.1)+
geom_smooth(data=quants,aes(x=rich.o, y=upper), se=F, size=0.1) +
geom_ribbon(data=quants,aes(x=rich.o, ymin=lower, ymax=upper,
y=mean), alpha=0.05, fill="blue",
size=0)+
geom_point(data=pontos, aes(x=richness.rarefaction, y=VPC, fill=abs(lat)),
col="black", size=2.5, pch=21)+
scale_fill_gradientn(colors=pal) +
scale_y_continuous(minor_breaks = NULL) +
scale_x_log10() +
geom_text(data=pvals, aes(x=x,y=y, label=sig), size=3, hjust=1,vjust=0,
col="black")+
xlab("Species richness (log10 scale)")+
ylab("VPC") +
theme_cowplot() +
theme(panel.background = element_rect(colour="lightgray"),
legend.position = "none")
fdiri_latpng(here("figures/FIG_S5.5_dirichlet_regressions_lat_exclude.png"), height = 600, width=800, res=100)
fdiri_lat
dev.off()## quartz_off_screen
## 2